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Abstract. We performed two-dimensional non-linear hydrodynamical simulations of two-component gravitating 
disks aimed at studying stability properties of these systems. In agreement with previous analytical and numerical 
simulations, we find that the cold gas component strongly affects the growth rates of the unstable global spiral 
modes in the disk. Already a five percent admixture of cold gas increases approximately two-fold the growth rate of 
the most unstable global mode while a twenty percent gas admixture enhances the modal growth rate by a factor 
of four. The local stability properties of a two-component disk coupled by self-gravity are governed by a stability 
criterion similar to Toomre's Q-parameter derived for one-component systems. Using numerical simulations, we 
analyse the applicability of a two-component local stability criterion for the analysis of the stability properties 
of global modes. The comparison of non-linear simulations with the analytical stability criterion shows that the 
two-component disks can be globally unstable while being stable to the local perturbations. The minimum value of 
the local stability criterion provides, however, a rough estimate of the global stability properties of two-component 
systems. We find that two-component systems with a content of cold gas of ten percent or less are globally stable 
if the minimum value of the stability parameter exceeds ~ 2.5. 

Key words. Galaxies: kinematics and dynamics - Galaxies: spiral 



1. Introduction 

Global spiral arms observed in disk galaxies are the re- 
sult of density waves propagating in the galactic disks 
(e.g. Binney & Tremaine 1987). Starting from pioneering 
works by B. Lindblad (1941, 1942), a number of authors 



extensively investigated this idea, and a variety of tech- 
niques have been developed to study the stability prop- 
erties of gravitating disks. Most of these studies concen- 
trated on one-component gaseous or coUisionless stellar 
disks. Galactic disks, however, are multi-component sys- 
tems whose properties depend on all the constituents. Lin 
& Shu (1966) first addressed the question showing the rela- 
tive importance of gas in providing the spiral gravitational 
field. Jog & Solomon ( 1984| ) (hereafter JS) considered the 
influence of a cold component on the stability properties of 
gravitating disks. They derived a local dispersion relation 
for a two-fluid instability of a star-gas disk and discussed 
a relative contribution of two components in the stability 
properties of the disk. In further investigations the results 
of JS were generahzed by Bertin & Romeo (1988) who 
studied the stability properties of the discrete global spiral 
modes by applying asymptotic WKB-methods. They also 
introduced the concept of the effective Q-parameter of the 
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two-component disks. Elmegreen (1995) re-considered the 
original equations derived by JS and reduced them to a 
single effective Q-parameter. Romeo (1992, 1994) included 
the effects of the finite thickness of stellar and gaseous 
disks. Recently, Rafikov (2001) gave the analysis of the 
axisymmetric gravitational stability of the thin rotating 
disk consisting of several components. 

The non-asymptotic global modal stability properties 
of multi-co mpon ent disks were studied to a lesser extent. 
Noh et al. ( 1991 ) studied the global stability properties of 
a thin two-fluid protoplanetary disk composed of gas and 
dust coupled by gravity and friction and found that the 
stability properties of these disks are strongly affected by a 
small amount of dust. Kikuchi et al. ( |1997| ) studied global 
stability properties of a three-phase disk with phases cou- 
pled by self-gravity and interchange processes between the 
phases. Korchagin & Theis (1999) numerically studied the 
behavior of global modes in multi-component disks at the 
non-linear stage of global instability. 

In this paper, we return to the problem of the global 
stability properties of multi-component disks. We consider 
the simplest possible system of that kind - a two-fluid thin 
disk coupled by the self-gravity of the components. We 
compare the results of our non-linear two-dimensional nu- 
merical simulations with the local stability criteria derived 
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by JS and study the applicability of the Jog-Solomon cri- 
terion for the prediction of the global stability properties 
of multi-phase gravitating disks. We also study the non- 
linear behaviour of the unstable global modes in multi- 
phase disks determining the dependence of the saturation 
levels of the unstable global modes on the mass fraction 
of a cold component. 

Sect. ^ describes the basic equations and the numer- 
ical code used in simulations. In Sect. ^ we present the 
equilibrium model used in our simulations. In Sect. ^ we 
briefly describe the two-fluid local stability criterion. Sect. 
^ presents the results of our simulations and gives a com- 
parison of the results of the numerical simulations with 
the two-fluid analytical stability criterion. 

2. Basic Equations 

To model the behaviour of the global perturbations in 
a gravitationally coupled multi-component disk we use 
the hydrodynamic approximation, modelling both a stel- 
lar and a gaseous component as polytropic fluids: 



(1) 



Here Ps^g are the vertically integrated pressures of each 
component, Kg^g are the polytropic constants and js,g are 
the polytropic indices for stellar and gaseous disks, respec- 
tively. In this approach, stellar and gaseous 'fluids' differ 
by the values of their polytropic indices js,g- We choose 
the value 7s = 2 for the stellar disk which reproduces 
the empirical "square root" proportionality between the 
velocity dispersion and the surface density found in galac- 



tic disks (Bottema 1992). Additional arguments justifying 



the description of the coUisionless stellar disks by a flui d 
dynamical approach are given by Kikuchi et al. ( 1997 ). 
The polytropic index for the gaseous disk is assumed to 
be 7g = 1.66. 

The behaviour of a two-fluid gravitationally interact- 
ing system is described by the continuity and momentum 
equations for each component: 



da, 



dt 



~{r'^.gUs,g) +]:^{.s,gVs^g) =0 



(2) 



dt ''^'^ dr r d(p r 



1 dPs^g d 



as,g dr dr 



(3) 



dVs,g ^ ^ dVs^g ^ Vs,g dVs,g ^ Vs,gUs,g 
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(4) 



Additionally, all components are coupled by their gravity, 
which can be expressed by the Poisson equation 



-Rout /•27T / / / ,/ 



{a,{r',^') + ag{r',cj,')ydr'd(t)' 
B.in Jo -y/r^ + r'2 — 2rr' cos{(j) — cj)') 



Here Ug g and v^ g are the radial and azimuthal velocities 
of stars and gas within the disk, and as^g is the surface 
density of each component. The potential derived from 
the self-gravity of a two-fluid disk, the explicit contribu- 
tion H from the rigid halo, and the pressure gradient 
determine the behaviour of the disk. All of the dependent 
variables are functions of the radial coordinate r, the az- 
imuthal angle 0, and the time t. The variables and the 
parameters of the model are expressed in units in which 
the gravitational constant G is equal to unity, the unit of 
mass is equal to 10^°Mq and the unit of length is 2 kpc. 
This choice determines the velocity unit Kinit = 149.1 
kms~^ and the time unit Tunit = 1-32 x 10^ yr. 

To solve the two-component hydrodynamical equa- 
tions (||) - we use the multi-component fluid- 
dynamical code developed by Korchagin & Theis (1999). 



.(5) 



The code is a multi-phase realization of the Eulerian 



ZEUS-2D code (Stone & Norman |1992| ) with the Van-Leer 
advection scheme. The two-fluid hydrodynamical equa- 
tions are solved using equally spaced azimuthal and log- 
arithmically spaced radial zones. In the numerical simu- 
lations discussed here we employed a grid of 256 x 256 
cells. The Poisson equation is solved by the 2D Fourier 
convolution theorem in polar coordinates. We look for the 
global modes starting from initial random perturbations 
at the level 10"^. 

3. Equilibrium model 

The equilibrium disk taken in our simulations is based 
on the equilibrium properties of real galaxies. Specifically, 
we applied the mo del o f the spiral galaxy NGC 1566 by 
Korchagin et al. (|2000D . The disk of NGC 1566 is ob- 
servationally well studied. Its rotation curve taken from 
Bottema ( 1992 ) can be approximated by 

Vir V2r 



Voir) 



(r2+i?2)3/4 (r2-^i?2)3/4 



(6) 



The radial dependence of the vertical velocity dispersion 
of the stars in the disk of NGC 1566 has an exponential 
distribution 



Cz = c^(r) = c^QCxp{-r/hcJ. 



(7) 



This uniquely determines the radial dependence of the ax- 
isymmetric surface density distribution assuming that the 
disk is self-gravitating in the vertical direction: 



cr(r) = cl{r)/TrGzo. 



(8) 



Here G is the gravitational constant, and zq is the effective 
thickness of the disk. According to equations (|^) and (||), 
the radial distribution of the stellar surface density can be 
expressed by 



(Ts{r) = (To^ exp(-r//i^)[l - (r/i?out) 



2i5 



(9) 



where is radial scale length for the stellar surface den- 
sity distribution. For numerical purposes we use in equa- 
tion (|^) a multiplier [1 — (r/_Rout)^]^ which results in a 



Orlova N., Korchagin V. and Theis Ch.: Global spiral modes in multi-component disks 



3 



Table 1. Parameters of the equilibrium model used in the 
numerical simulations 



Inner disk radius (Rin), kpc 


0.2 


Outer disk radius (i?out), kpc 


10.0 


Disk thickness (zo), kpc 


0.7 


Scale length of velocity dispersion {hc^), kpc 


2.6 


Scale length of surface density (ha), kpc 


1.3 


Central velocity dispersion (c^q), kms~^ 


168 


Total mass of the disk (Md), Mq 


1.78 X 10"' 


Rotation curve parameters: 




Vi, V2, kms-i 


399, 82 


Ri, R2, kpc 


4.4, 0.6 



vanishing surface density at the outer boundary of the 
disk. We further assume that the surface density profile of 
the gaseous disk and its thickness are the same as those 
in the stellar disk. However, the total mass of the gaseous 
disk is smaller compared to the mass of the stellar disk 
and varies in our models. 

The radial velocity dispersion of the stellar and gaseous 
components can be derived from the equation of state by 



-Ul/2 



(10) 



The free parameters i^s.g can be fixed by the equations 
(0) and and by the assumption of a constant ratio of 
the vertical and radial velocity dispersions. For the lat- 
ter, to get the disk with weak global modes we have cho- 
sen Csz/csr = 0.9. Our choice can be justified as follows. 
According to Bottema (1992), the measured velocity dis- 
persion in the disk of NGC 1566 gives the real dispersion 
perpendicular to the galactic disk. In our one-component 
model, we fix the ^-component of the velocity disper- 
sion and hence the surface density and the total mass 
of the disk, and allow for a variation of the radial ve- 
locity dispersion. Increasing the ratio Csz/Cgr we decrease 
the radial velocity dispersion, hence destabilizing the disk. 
Simulations show that for the measured z-dispersion and 
rotation curve in the disk of NGC 1566, the disk is stable 
for the usually assumed ratio Cgz/csr = 0.6. Therefore, we 
choose in our studies the ratio Csz/csr = 0.9 which leads to 
a weakly unstable disk. This ratio is larger than the mea- 
sured ratio of the velocity dispersions Csz/csr = 0.5 — 0.6 in 
solar neighborhood, and close to the upper limit measured 
in other galaxies. E.g., the determination of this ratio in 
the disk of galaxy NGC 488 undertaken by Gerssen et al. 
(1997) gives the value 0.7±0.19. With our choice, however, 
we do not aim to model any particular galaxy. 

The equilibrium rotation of the disk given by the equa- 
tion (^) is balanced by the self-gravity, the pressure gradi- 
ent of the disk and the external potential of a rigid halo: 



r a dr dr \ 



H 



(11) 



By this, we can fix the halo properties. Knowing 
the initial surface density of the disk we get the disk's 
self-gravity by solving the Poisson equation numerically. 
Equation ( [ll] ) can then be used to calculate the gradient 



of the halo potential and the halo mass distribution. The 
general parameters of our equilibrium model are listed in 
Table 0. 



4. Stability Criteria 

The local stability of a thin gravitating disk can be ex- 
pressed in terms of Toomre's Q parameter which for the 
gaseous disk has the form (Safronov |196|, Toomre p^e^ : 



Q 



c{r)K{r) 
1:0(7 [r) 



(12) 



Here, c(r) is the sound speed of the gas, K{r) is the 
epicyclic frequency and u{r) is the surface density of the 
disk. A gravitating disk is locally unstable when Q < 1. 

The stability criterion (|l2|) was derived for tightly 
wound perturbations. However, it can be successfully used 
as an indicator for the stability of global modes in grav- 
itating disks. Numerical simulations show that a locally 
stable stellar or gaseous disk can be globally unstable if 
the minimum value of Q does not exceed ^ 2. For exam- 
ple, a one-component disk with the rotation and surface 
density distribution given by equations (^) and (^) can 
be globally unstable, if the minimum value of Q does not 
exceed ~ 1.6 (Korchagin et al. 2000| ). 

Disks of real galaxies consist of stellar and gaseous 
components. The stellar component comprises the largest 
fraction of the total mass of the disk. However, colder 
gaseous components can play a significant role for the sta- 
bility properties of multi-component disks. JS derived a lo- 
cal stability criterion for a two-component fluid disk which 
can be written in the form Qjs < 1 where 



gjs = 



(13) 



Here Qs and Qg are Toomre's Q-parameters deter- 
mined for the 'stellar' and gaseous fluids, correspondingly. 
q is the dimensionless radial wavenumber of perturbations 
q = kcs/n. The relative contribution F of gas to stars is 
given by the ratio of the second to the first term of the 
stability parameter (113): 



F = 



1 



1 + q^cl/cl 



(14) 



It can be seen from equation (^3|) that the relative 
contribution of the components is always larger than the 
ratio of the surface densities of the components, provided 
the gas is dynamically colder than the stars. In the limiting 
case of large wavenumbers g — > cxa the contribution of the 
components is given by the ratio UgC^/a^c^ - the result. 



obtained by Marochnik & Suchkov (1974). 

We will use the criterion given by equation ( [l3| ) for 
the analysis of instability of the global modes in two- 
component gravitating disks. 
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Fig. 1. Equilibrium curves for the one-component 'stel- 
lar' disk. Angular velocity to and epicyclic frequency k are 



given in units of 0.75 x 10~ yr~ . 
is given in units 2.5 x IO^MqIcpc" 
length L = 2 kpc. 



5. Results 



The surface density a 
^, radius is in units of 



To study the influence of a cold component on the overall 
dynamics of a two-fluid gravitating disk, we start our anal- 
ysis with a one-component 'stellar' disk which we choose 
to be close to marginal stability. Fig. ^ shows the equilib- 
rium properties of this disk. The local stability parameter 
Q is well above unity with a minimum value Qmin ~ 1-6. 

The radial behavior of the Q-parameter is one of the 
key assumptions in the modal analyses of spiral density 
waves. In previous studies it was usually assumed that 
the profile of the stability parameter Q has a barrier at 
the central regions of the disk and is equal to unity in the 
outer regions of the disk (e.g. Bertin et al. 1989). Such 
a behavior of an 'effective' Q-parameter is supposed to 
model the absence of cold gas as well as the transition in 
geometry from the disk to the nuclear bulge in the central 
regions of a galaxy, and to describe the increasing role of 
a cold gaseous component in the outer regions of a galac- 
tic disk. In this paper we do not introduce an effective 
Q-parameter, but we try to base our studies on obser- 
vational properties of disk galaxies. We choose the basic 
state in our fiducial model based on the properties of the 
nearby late type spiral galaxy NGC 1566. The 'stellar' 
Q-parameter shown in Fig. |l| has an outer barrier which 
is typical for disks with exponentially decreasing radial 
surface density distributions and exponentially decreasing 
velocity dispersions. 

Linear global modal analysis of this model reveals a 
few slowly growing modes with the fastest growing spi- 
ral to — 2. Nonlinear 2D simulations confirm this result 
giving an exponential growth rate of the m = 2 mode of 
approximately 0.11 in dimensionless units. Linear theory 
predicts a value of 0.12 which is in good agreement with 
non-linear simulations. Fig. illustrates the evolution of 




Fig. 2. Temporal evolution of the to — 1—4 global Fourier 
amplitudes computed for the one-component 'stellar' disk 
of Fig. H Time is in units of 1.32 x 10'' yr. 

the to = 1 — 4 global modes in the one-component 'stellar' 
disk given in terms of the global Fourier amplitudes 



1 



2ir 



(15) 



The most unstable to = 2 mode grows linearly until 
t « 110, and then saturates at a level A2 « 0.2 (Fig. ||). 
The linear regime corresponds to an exponential growth 
of perturbations. 

In the next step we split the one-component disk into a 
'stellar' and a 'gaseous' subsystem keeping the total mass 
of the disk constant. We consider a set of models with a 
cold component consisting of 5%, 10% and 20% of the to- 
tal mass of the disk. In the first series of simulations we 
keep the 'stellar' and gaseous equations of state the same 
for all the models. As a result, the gas velocity dispersion, 
which depends on the gas fraction, is not constant. Thus, 
the velocity dispersion of gas is about 20% of the velocity 
dispersion of the stars in the model with 5% gas admix- 
ture. This ratio increases up to 36% for the model with 
20% gas admixture. 

The strong influence of the gas fraction is shown in 
Fig. 3 which displays the growth rate of the most unstable 
TO = 2 mode in a one-component 'stellar' disk compared 
to the modal growth rate of the two-component models 
with 5%, 10%, and 20% gas admixtures. Already a five 
percent gas fraction enhances the modal growth rate by a 
factor of 1.7 compared to the one-component disk. A 20% 
gas admixture results in an increase of the growth rate by 
a factor of 3.7. Table § gives the masses and the dispersion 
ratio of the components in the computed models together 
with the growth rates of the most unstable modes, satura- 
tion levels and the minimum values of the two-component 
stability criterion Qjg calculated by equation ([l^). The 
growth rate of the dominant mode increases in the mod- 
els with higher gas content, being the same for the gaseous 
and the stellar components. However, the saturation level 
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Fig. 3. Time dependence of the global Fourier amplitudes 
for the most unstable m = 2 spiral mode computed for 
models with different gas fractions (I: 5%, II: 10%, III: 
20%). Time is in units of 1.32 x lO'' yr. 

Table 2. 











Saturation level 




Model 


MJMs 


Cg/Cs 


Im2 


gas 


stars 


QjSmin 


I 


0.05 


0.2 


0.18 


0.35 


0.27 


2.2 


II 


0.1 


0.25 


0.24 


0.32 


0.25 


1.75 


III 


0.2 


0.36 


0.4 


0.29 


0.22 


1.36 



Note. Dependence of the growth rates and the saturation 
levels on the content of a cold component for the governing 
m — 2 mode. The total mass of the disk is hxed to 1.78 x 
IO^^A/q. The growth rates are given in units of 0.75 x 

10-V"'- 



of the most unstable mode is systematically higher for 
the gaseous disk. In general, the saturation level slightly 
decreases in the models with higher gas content. A possi- 
ble explanation of this behavior is that spiral shocks play 
a more important role in the saturation of exponentially 
growing modes when more cold gas is present. 

The ratio of the velocity dispersions of the subsystems 
can also affect the stability of the disk. The influence of 
the velocity dispersions on the growth rates of the global 
modes was studied in a series of numerical experiments 
in which the fraction of the cold component was fixed 
to 10% of the total mass of the disk, while the velocity 
dispersion ratio was allowed to change. The results of these 
simulations are shown in Fig. 4. The growth rate of the 
most unstable global mode is only slightly affected by the 
ratio of the velocity dispersions. The saturation amplitude 
remains approximately at the same level. Therefore, the 
mass fraction of the cold component is a more important 
parameter for the growth rate of perturbations than the 
ratio of the velocity dispersions. This is in agreement with 
results of Bertin & Romeo (1988). They demonstrated (see 
their Fig. 5) that the criterion of marginal stability of a 
two-component disk also depends slowly on the ratio of the 



Fig. 4. Time evolution of the global Fourier amplitude 
for the most unstable mode (m = 2) computed for models 
with 10% gas admixture and a ratio Crs/crg of the velocity 
dispersions equal to 1.95, 3.9 and 7.8. Time is in units of 
1.32 X 10^ yr. 

velocity dispersions of the components and more strongly 
on the ratio of their densities. 

Disks with gaseous admixture show qualitatively the 
same evolution as the pure stellar disk. However, they 
evolve on a shorter time scale. The appearance of the 
global modes in all models is quite similar despite the 
considerable difference in their growth rates. Fig. ^ shows 
the logarithmically spaced surface density contour plots 
taken at the moments when the global amplitude of the 
TO = 2 mode in all four models is A2 « 0.06. All the models 
develop an open two-armed spiral of comparable winding. 
Note, however, that the spiral arms are more narrow, more 
tightly wound and less extended when more gas is present. 
Analysis of the global modes with help of the local disper- 
sion relation (Berlin & Romeo 1988) also shows a shift 
towards tighter spirals with a smaller corotation radius in 
models with higher gas content. 

To use the stability criterion ([Tsl), one needs to know 
the wavenumber of the perturbations. Fig. ^ shows the ra- 
dial dependence of the radial wavenumber for m = 2 per- 
turbations computed for the models I - III. The wavenum- 
bers were determined numerically by calculating the pitch 
angles of the spiral arms shown in Fig. |^ at different radii, 
i.e. we used the relation 



cot i = 



kR 

TO 



(16) 



To find a pitch angle of the spiral at a particular radius, 
we determine the line of the maximum amplitudes of the 
spiral arm and calculate the angle between the tangents 
to the line of maxima and the corresponding circle. 

Fig. shows the radial profile of the dimensionless 
Jog-Solomon stability parameter Qjs computed with the 
wavenumbers taken from Fig. |[ Fig. |^ further illustrates 
the destabilizing role of a cold gas component: the mini- 
mum value Qjs of the Jog-Solomon parameter decreases 
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-1 




-3-2-11 2 3 

Model III. Time=25 



-3 



-2 



-1 



Fig. 5. Surface density perturbations of the stellar component. Displayed are the spiral perturbations in the one- 
component model (upper left), and in the two-component models with different gas fractions. The snapshots are taken 
at the moment when the m = 2 global amplitude reaches A2 « 0.06. The contour levels are logarithmically spaced 
between the maximum value of the density perturbation and one-hundredth of the maximum value. Circles show the 
positions of the outer Lindblad resonance (dashed line) and corotation (solid line). 
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Fig. 6. The radial dependence of the radial wavenumber 



Fig. 7. The radial dependence of the local two-fluid sta- 

different gas 



of two-armed spiral for the models I - III. Unit of length ^ility parameter Qjg for disk models with 
L = 2 kpc fraction. Unit of length L = 2 kpc. 



with increasing content of cold gas. The minimum value minimum value of 2.2 found for the two-component model 

of the Jog-Solomon parameter Qjs is about 1.36 for the with 5% gas fraction. In all three models, the Jog-Solomon 
two-component system with 20% gas admixture versus a stability criterion exceeds unity, and all the models are 
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Fig. 8. The growth rate as a function of the minimum 
value of local stability parameter Qjs. 



locally stable. Such a behaviour is quite similar to the be- 
haviour of one-component disks which are known can be 
globally unstable even if the minimum value of Toomre's 
Q parameter is above unity. 

Decreasing the total mass of the disk results in its sta- 
bilization which consequently leads to an overall growth of 
the Jog-Solomon parameter Qjs- Fig. ^illustrates this by 
displaying the growth rate of the dominant mode to = 2 
as a function of the minimum value of the Jog-Solomon 
parameter Qjs- The latter was obtained from numerical 
simulations of two-component disks with 10% gas admix- 
ture. These disks are stabilized when the minimum value 
of the Jog-Solomon parameter exceeds ~ 2.5. 



6. Summary 

In this paper we studied the applicability of the two-fluid 



local stability criterion derived by Jog & Solomon (1984) 
for the onset and growth of non-axisymmetric perturba- 
tions in two-component gravitating disks with an open 
spiral structure. We performed three sets of numerical ex- 
periments: one set with fixed total mass of the disk and 
fixed polytropic equation of state, but variable gas frac- 
tion. Another set of simulations was performed with a 
fixed total mass of the disk and a fixed ratio of the stellar 
and gaseous velocity dispersions. A final set of calculations 
deals with a variable disk mass. Our numerical experi- 
ments have shown that the mass contribution of the cold 
component strongly influences the stability properties of 
the disk. Even a low gas fraction of 5% destabilizes the 
disk substantially, e.g. the growth rate of m = 2 pertur- 
bations increases by a factor 1.7. All our disk models are 
unstable to m = 2 spiral modes. More unstable disks de- 
velop narrower and more tightly wound spiral arms. The 
saturation levels are almost constant for all simulations, 
with a slightly enhanced level for the gaseous component. 
A comparison of the results of the direct two-dimensional 
numerical simulations with the JS stability criterion allows 



us to conclude that similarly to Toomre's stability param- 
eter Q derived for one-component disks, the Jog-Solomon 
stability parameter Qjs can be used for the prediction 
of the global stability properties of two-component grav- 
itating disks which greatly simplifies the analysis of disk 
stability. 
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